%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Empirical Application(Timber Dataset)                                                     %
% Test: Choice of Specification for the Bidders' Private Values (Linear versus Exponential) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all

load('data.mat'); % load timber dataset
vol2=vol2./1000;
vol3=vol3./1000;
apv2=apv2./10;
apv3=apv3./10;

L2=size(wb2,1); % number of auctions with N=2
L3=size(wb3,1); % number of auctions with N=3
L=L2+L3;

% covariates
x2=[ones(L2,1),apv2,vol2];
x3=[ones(L3,1),apv3,vol3];
x = [x2;x3];
p=size(x,2);

% winning bids
wb=[wb2;wb3];

N=[2,3]; % Number of bidders
alpha=0.12:0.02:0.80; % quantile levels alpha
psi2=ones(L2,1)*(N(1)*(alpha.^(N(1)-1))-(N(1)-1)*(alpha.^N(1))); % Psi function given N=2
psi3=ones(L3,1)*(N(2)*(alpha.^(N(2)-1))-(N(2)-1)*(alpha.^N(2))); % Psi function given N=3
psi = [psi2;psi3];
m=size(alpha,2);

B=5000; % Bootstrap replications
toler=0.1; % Parameter of tolerance for precision of the estimation

gammai=[1,0.3,0.3]'; % Initial guess for the parameters (chosen to speed up the estimation)

dat=[wb,x]; % Matrix with data

    % Computation of the choice of specification test statistic 
    for k=1:m  
       % Obj F under H0
       [objE,~,~]=MM_Exp(psi(:,k), p, wb, x, toler, gammai);
       % Obj F under H1
       [objL,~,~]=MM(psi(:,k), p, wb, x, toler, gammai);
       Vuong(k,:)=objL-objE;
    end

    Vuongmax=sqrt(L)*max(Vuong);  % test statistic
     
% Pairwise Bootstrap of the choice of specification test statistic
for b=1:B
    datb2 = dat(randsample(1:L2,L2,true),:); % Resampling of the matrix of data in the subsample for N=2
    datb3 = dat(randsample(L2+1:L,L3,true),:); % Resampling of the matrix of data in the subsample for N=3
    wb2 = datb2(:,1); % Resampled winning bids N=2
    wb3 = datb3(:,1); % Resampled winning bids N=3
    wb=[wb2;wb3];
    x2 = datb2(:,[2 3 4]);  % Resampled covariate N=2
    x3 = datb3(:,[2 3 4]); % Resampled covariate N=3
    x=[x2;x3];
    
    for k=1:m
        % Obj F under H0
        [objEb,~,~]=MM_Exp(psi(:,k), p, wb, x, toler, gammai);
        % Obj F under H1
        [objLb,~,~]=MM(psi(:,k), p, wb, x, toler, gammai);
        Bootvuong(k,:)=sqrt(L)*((objLb-objEb)-Vuong(k,:));
    end
   
   Bootvuongmax(:,b)=max(Bootvuong);
   b        
end

% Critical Values
c99=quantile(Bootvuongmax',0.99);
c95=quantile(Bootvuongmax',0.95);
c90=quantile(Bootvuongmax',0.90);

% P-value computation    
R=zeros(1,B);
for b=1:B
    if Bootvuongmax(:,b)>= Vuongmax
            R(:,b)=1;
    else
            R(:,b)=0;
    end
end
    pvalue=sum(R')/B